---
title: "Individual-Level Analysis"
author: "Sally Sharif"
date: "2026-01-18"
output:
  word_document: default
  html_document: default
---

```{r}
rm(list = ls())

# loading required packages
library(table1)
library(corrtable)
library(dplyr)
library(lmtest)
library(texreg)
library(nnet)  
library(sandwich)
library(clubSandwich) 
library(janitor)  
library(VGAM)

# loading the data 
obj_name <- load("individual_level_data.RData")
data_ind <- get(obj_name)
```


```{r descriptive stats}

# Appendix Table 11: descriptive statistics 
table1(~ factor(report_state) + factor(report_military_police) + factor(report_local) +
         factor(protection) + factor(collective) + factor(community) + factor(govsupport) +
         factor(internationalsupport) + age + factor(gender) + 
         factor(education) + factor(friends) + factor(married) + 
         factor(employed) + factor(disability), data=data_ind)

# Appendix Table 12: correlation table
corrdata <- data_ind %>% 
  dplyr::select(report_state, report_military_police, report_local, protection, 
                collective, community, govsupport, internationalsupport, 
                employed, age, gender, education, friends, married, disability) 

corrdata_ok <- corrdata %>%
  select(where(~ all(is.na(.x)) == FALSE)) %>%
  select(where(~ sd(.x, na.rm = TRUE) > 0))

nrow(corrdata_ok) # 4,605 because all ex-combatants that are not in productive projects are removed

correlation_matrix(corrdata_ok, digits = 2, use = "upper", replace_diagonal = TRUE)
```

```{r Table 3}
## Table 3: reporting security threats to the state (columns 1-2) 

# (1) State, simple
m1_state_simple <- glm(report_state ~ protection,
                       data = data_ind, family = binomial())
coeftest(m1_state_simple, vcov = vcovHC(m1_state_simple, type = "HC1"))

# (2) State, with controls
m2_state_full <- glm(report_state ~ protection + employed + age + gender + 
                       education + friends + married + disability,
                       data = data_ind, family = binomial())
coeftest(m2_state_full, vcov = vcovHC(m2_state_full, type = "HC1"))

# Province-clustered SEs
dat_cl <- subset(data_ind, !is.na(departamento))
m2_state_full_cl <- glm(report_state ~ protection + employed + age + gender + 
                          education + friends + married + disability,
                          data = dat_cl, family = binomial())
coeftest(m2_state_full_cl, vcov = vcovCL(m2_state_full_cl, 
                                              cluster = ~ departamento, type = "HC1"))

## Table 3: reporting security threats to military/police (columns 3-4) 

# (3) Military/police, simple
m3_mil_simple <- glm(report_military_police ~ protection,
                     data = data_ind, family = binomial())
coeftest(m3_mil_simple, vcov = vcovHC(m3_mil_simple, type = "HC1"))

# (4) Military/police, with controls
m4_mil_full <- glm(report_military_police ~ protection + employed + age + gender +
                     education + friends + married + disability,
                   data = data_ind, family = binomial())
coeftest(m4_mil_full, vcov = vcovHC(m4_mil_full, type = "HC1"))

# Province-clustered SEs (drop NA clusters)
m4_mil_full_cl <- glm(report_military_police ~ protection + employed + age + gender +
                        education + friends + married + disability,
                      data = dat_cl, family = binomial())
m4_mil_CL <- coeftest(m4_mil_full_cl,
  vcov = vcovCL(m4_mil_full_cl, cluster = ~ departamento, type = "HC1"))

# Table 3: reporting security threats to non-state actors (columns 5-6) 

# (5) Non-state, simple
m5_loc_simple <- glm(report_local ~ protection,
                     data = data_ind, family = binomial())
coeftest(m5_loc_simple, vcov = vcovHC(m5_loc_simple, type = "HC1"))

# (6) Non-state, with controls
m6_loc_full <- glm(report_local ~ protection + employed + age + gender +
                     education + friends + married + disability,
                   data = data_ind, family = binomial())
coeftest(m6_loc_full, vcov = vcovHC(m6_loc_full, type = "HC1"))

# Province-clustered SEs (drop NA clusters)
m6_loc_full_cl <- glm(report_local ~ protection + employed + age + gender +
                        education + friends + married + disability,
                      data = dat_cl, family = binomial())
coeftest(m6_loc_full_cl, vcov = vcovCL(m6_loc_full_cl, cluster = ~ departamento, type = "HC1"))

# Table 3
texreg(list(m1_state_simple, m2_state_full, m3_mil_simple, m4_mil_full, m5_loc_simple, m6_loc_full))
```

```{r Table 4}
### Reporting behavior when in collective projects or collaborating with civilians
## Table 4: reporting security threats to the state (columns 1-2) 

# (1) State ~ Collective + Collective×Protection + controls
m1_state_collect <- glm(
  report_state ~ collective + protection + collective:protection +
    age + gender + education + friends + married + disability,
  data = data_ind, family = binomial())

coeftest(m1_state_collect, vcov = vcovHC(m1_state_collect, type = "HC1"))

# with province clusters
m1_state_collect_CL  <- glm(
  report_state ~ collective + protection + collective:protection +
    age + gender + education + friends + married + disability,
  data = dat_cl, family = binomial())

coeftest(m1_state_collect_CL, vcov = vcovCL(m1_state_collect_CL, cluster = ~ departamento, type = "HC1"))

# (2) State ~ Community + Community×Protection + controls
m2_state_comm <- glm(
  report_state ~ community + protection + community:protection +
    age + gender + education + friends + married + disability,
  data = data_ind, family = binomial())

coeftest(m2_state_comm, vcov = vcovHC(m2_state_comm, type = "HC1"))

# with province clustered SEs
m2_state_comm_CL <- glm(
  report_state ~ community + protection + community:protection +
    age + gender + education + friends + married + disability,
  data = dat_cl, family = binomial())

coeftest(m2_state_comm_CL, vcov = vcovCL(m2_state_comm_CL, cluster = ~ departamento, type = "HC1"))

# Reporting to military/police (Columns 3–4)

# (3) Mil/Police ~ Collective + Collective×Protection + controls
m3_mil_collect <- glm(
  report_military_police ~ collective + protection + collective:protection +
    age + gender + education + friends + married + disability,
  data = data_ind, family = binomial())

coeftest(m3_mil_collect, vcov = vcovHC(m3_mil_collect, type = "HC1"))
 
# province clusters
m3_mil_collect_CL <- glm(
  report_military_police ~ collective + protection + collective:protection +
    age + gender + education + friends + married + disability,
  data = dat_cl, family = binomial())

coeftest(m3_mil_collect_CL, vcov = vcovCL(m3_mil_collect_CL, cluster = ~ departamento, type = "HC1"))

# (4) Mil/Police ~ Community + Community×Protection + controls
m4_mil_comm <- glm(
  report_military_police ~ community + protection + community:protection +
    age + gender + education + friends + married + disability,
  data = data_ind, family = binomial())

coeftest(m4_mil_comm, vcov = vcovHC(m4_mil_comm, type = "HC1"))

# with province clustered SEs
m4_mil_comm_CL <- glm(
  report_military_police ~ community + protection + community:protection +
    age + gender + education + friends + married + disability,
  data = dat_cl, family = binomial())
  
coeftest(m4_mil_comm_CL, vcov = vcovCL(m4_mil_comm_CL, cluster = ~ departamento, type = "HC1"))

# Reporting to non-state actors (Columns 5–6)

# (5) Non-State ~ Collective + Collective×Protection + controls
m5_local_collect <- glm(
  report_local ~ collective + protection + collective:protection +
    age + gender + education + friends + married + disability,
  data = data_ind, family = binomial())

coeftest(m5_local_collect, vcov = vcovHC(m5_local_collect, type = "HC1"))

# clusters
m5_local_collect_CL <- glm(
  report_local ~ collective + protection + collective:protection +
    age + gender + education + friends + married + disability,
  data = dat_cl, family = binomial())
  
coeftest(m5_local_collect_CL, vcov = vcovCL(m5_local_collect_CL, cluster = ~ departamento, type = "HC1"))

# (6) Non-State ~ Community + Community×Protection + controls
m6_local_comm <- glm(
  report_local ~ community + protection + community:protection +
    age + gender + education + friends + married + disability,
  data = data_ind, family = binomial())

coeftest(m6_local_comm, vcov = vcovHC(m6_local_comm, type = "HC1"))

# province clusters
m6_local_comm_CL <- glm(
  report_local ~ community + protection + community:protection +
    age + gender + education + friends + married + disability,
  data = dat_cl, family = binomial())

coeftest(m6_local_comm_CL, vcov = vcovCL(m6_local_comm_CL, cluster = ~ departamento, type = "HC1"))

# Table 4
texreg(list(m1_state_collect, m2_state_comm, m3_mil_collect, m4_mil_comm, m4_mil_comm_CL, m5_local_collect,m6_local_comm))
```

```{r Appendix Table 23}
## Multinomial logit models of reporting behavior

# building overlap indicators (state_any vs nonstate_any)
# treating "state" as (report_state == 1 OR report_military_police == 1)
# treating "non-state" as (report_local == 1)
dat <- data_ind %>%
  mutate(
    state_any = as.integer(report_state == 1),
    nonstate_any = as.integer(report_local == 1),
    both_state_nonstate = as.integer(state_any == 1 & nonstate_any == 1),
    neither = as.integer(state_any == 0 & nonstate_any == 0),
    cat4 = case_when(
      state_any == 1 & nonstate_any == 0 ~ "State only",
      state_any == 0 & nonstate_any == 1 ~ "Non-state only",
      state_any == 1 & nonstate_any == 1 ~ "Both",
      TRUE                                ~ "Neither")) %>%  mutate(
    cat4 = factor(cat4, levels = c("Neither","State only","Non-state only","Both")))

# Appendix Table 22: Descriptive share of overlap 
overlap_tab <- dat %>%
  count(state_any, nonstate_any) %>%
  mutate(share = 100 * n / sum(n))

overlap_tab

# fitting multinomial model (base category = "Neither")
mnl_vglm <- vglm(
  cat4 ~ protection + employed + age + gender + education + friends + married + disability,
  family = multinomial(refLevel = "Neither"),
  data = dat)

# Appendix Table 23
texreg(mnl_vglm)
```

